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Abstract 



In the framework of NRQCD and parton model, we estimate in detail, the production 
cross section for polarized J/ip from B meson decays. In order to contrast with data, 
we also take into account additional J/ip production due to decay of excited charmonium 
states. We calculate the helicity parameter, a, and as an application, we study our results 
for the Tevatron. This is in contrast to the earlier studies which were performed for 
prompt J/ip production from pp collisions. Our estimates are, for J ftp from B decays, 
c^j/tp = —0.04 ± 0.06 and for B decays to ip', c^V' ~ —0.03 ± 0.07. These results have been 
evaluated in the J/ip transverse momentum interval, 10 GeV < /ct < 30 GeV. In the limit 
of the color singlet model, a shows a direct dependence on the Peterson parameter, thereby 
reflecting the dynamics of the b quark hadronization. With Run II of the Tevatron, it is 
expected that the fits for a will improve by about a factor of 50, leading to better limits 
on the matrix elements. 
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1 Introduction 



The production mechanism of bound states involving a heavy quark and anti-quark system 
can be addressed within non-relativistic quantum chromodynamics (NRQCD) In the 
earhest attempts, charmonium production was described by the color singlet model through 
processes like B decays (6 ccs) and gluon-gluon fusion {gg — > ccg) p. We refer 

to for a review on these issues. However, the color singlet models had several problems, 
e.g., underestimation of the hadropro duct ion of charmonium 0, the ip' anomaly P,^] and 
infrared divergences in P wave charmonium production, which later had a resolution based 
on factorization results . These problems suggested the need to advance beyond the 

color singlet model or similar variants like the color evaporation model [0. In a systematic 
approach, by including the color octet contributions within the NRQCD framework it was 
shown that these problems could indeed be resolved to a good accuracy . 

Within NRQCD, which is well designed for separating relativistic from non-relativistic 
scales, Bodwin, Braaten and Lepage developed a factorization formalism to calculate 
quarkonium decays and production |jlO[ . The formalism allows for a systematic calculation 
of the inclusive cross sections to any order in strong coupling a^, and an expansion in v"^. 
Here, v is the relative velocity of the quark and antiquark and is inversely proportional to 
the heavy quark mass. As an illustration, following potential model calculations, for botto- 
nium states, f ~ 0.1 while for charmonium states, v ~ 0.3 indicating a better convergence 
in the perturbative expansion for heavier quark states [jl6|. It is interesting to observe 



that the theory exhibits a scale hierarchy of the type, mq 3> mqv ^ mqv ~ Aqcd, where 
mq is the heavy quark mass. Therefore, it is appropriate to use NRQCD as an effective 
field theory with v as the expansion parameter which is also a naturally small scale of the 
theory. In addition, at leading order in v, NRQCD has a strong correspondence to a 1/mq 
expansion as in heavy quark effective theory. 

A large fraction of the existing literature on NRQCD employing the factorization for- 
malism, broadly concentrates on one of the following two issues; (i) on prompt charmonium 
production in hadron |T^,|T^, 7p and ep and in e'^e" 

or (ii) on charmonium production in hadronic B decays [p5|,p^. 



^0^0 collisions 
Prompt production 

refers to quarkonium (here charmonium) that is created in interactions of the colliding 
particles or their constituents, while charmonium production is also possible in weak de- 
cays of B mesons, which will be the focus of the present work. At the functional level, the 
calculations in case (i) have been adopted to phenomenologically extract NRQCD matrix 
elements from experimental data on charmonium production. This is possible, because 
these calculations incorporate bound state effects of the initial hadronic states, usually in 
the framework of the QCD improved parton model (PM). In the case of semi-inclusive B 
decays with charmonium final states, the ACCMM model |]28| and the PM have been 
successfully adopted for this description. A central feature of these results has been to 
illustrate the importance of color octet elements to accommodate the observed momentum 
spectra of J/i^. 

Quarkonium polarization provides an additional test of the color octet production mech- 
anism of NRQCD pOi. The polarized cross section has been calculated for prompt ip pro- 
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duction [0 as well as for J/ip production in h quark decays We remark that, in 

the case of J/ ip production from h decays, the calculations do not take into account bound 
state effects, although they have been estimated to be significant [32]. Hence, a comparison 
with data in this case may not be too meaningful, given the uncertainties, originating from 
the negligence of the initial hadron, and the additional errors due to the non-perturbative 
NRQCD matrix elements. As a salient prediction, within NRQCD, prompt charmonium 
production is expected to be predominantly in transverse polarization state for large trans- 
verse momenta {px) |13, 1^,18,^; but this prediction is not in agreement with the CDF 



data ||3J]. We note that the polarization prediction arises from the dynamics of massless 
partons and for large pr, the role of gluon dynamics is important in prompt quarkonium 
production, especially through the dominance of gluon fragmentation Furthermore, 
a gluon couples easily to the ^Si color octet state, which is expected to be a dominant 
spectral state in the prompt J/ ip production mechanism at the Tevatron. But correspond- 
ingly, the charmonium production at large transverse momenta (with > 20 GeV) are 
not fully probed by current experiments. Besides, there are large errors in the polarization 
measurements. Therefore, these features alone preclude any possible conclusions on the 
predictions by NRQCD for polarized charmonium production. 

On the other hand, at moderate transverse momenta (with p^ < 20 GeV), one can 
perform the polarization studies for prompt charmonium production to make an estimate 
of the color octet elements and also compare with the unpolarized cross sections. This is of 
particular relevance to the Tevatron where there are no complications due to higher twist 
effects ||3^. We refer to ||3^ for an update on prompt production of polarized charmonium 
for the Tevatron. Simultaneously, one can also study the charmonium production which is 
not prompt and in particular estimate the cross sections for polarized production. 

It is the goal of the present work to analyze the polarization predictions for J/ip from 
B meson decays at the Tevatron. Unlike in the case of prompt production, in this process, 
we do not expect gluon fragmentation as a dominant source for J/ip production, which 
led to predominantly transverse polarized J/ip. Therefore, our calculation can serve as an 
independent probe of NRQCD dynamics for polarized charmonium production, besides the 
existing knowledge from prompt production. We employ the PM approach as discussed 
p9|] to fold the quark level calculations to arrive at a i? hadron decay. We calculate the 



m 



helicity parameter, a, from the production cross section of the three polarization states 
of the J/V'- We observe that a significant drawback of any such analysis are due to our 
present poor understanding of the relevant NRQCD matrix elements. As an outcome of our 
approach, we note that in the color singlet model (when the octet elements are set to zero), 
the a prediction reduces to the details of bound state effects of the PM. In other words, 
a depicts a strong dependence on the Peterson fragmentation function which describes 
the Fermi motion of the h quark in the B meson. In some sense, this result is also to be 
anticipated simply on grounds that the color singlet model predictions depend on the shape 
and form of the initial state wave function of the decaying system. With data from Run 
II of Tevatron, which is expected to increase the accuracy by a factor of 50, our analysis 
may be useful to tighten the estimates for the matrix elements significantly and make 
our polarization estimates more precise In addition, in the future, a complete global 
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fit/analysis to quarkonium production will certainly make the predictions more robust. 

Our paper is organized as follows. In the next section, for completeness, we review 
the basic ideas of NRQCD pertinent to our calculations. In section ^ we introduce the 
effective Hamiltonian which describes the quarkonium production process through free b 
quark decays. Using this formalism, we study the semi-inclusive decay of a free b quark 
into J/ip. The short-distance coefficients and the NRQCD matrix elements are explicitly 
calculated and the decay width is presented for b J/i/j{X) + X. Here, A denotes one of the 
three helicity states of J/ip. Towards the end of this section, we also discuss an extension 
of our calculations to excited charmonium states. The bound state effects of the initial 
B meson, whose influence hitherto has been neglected in the calculation, are described in 
section |[ In this analysis, we use the PM approach to evaluate the bound state effects. 
Starting with a short introduction to the PM, the restrictions of the model applicability 
and estimates for the semi-inclusive decay rate for a B meson with a charmonium flnal 
state are presented. This is followed by section ^ where we describe the application of 
the results derived so far to the Tevatron and introduce suitable kinematic variables. In 
section |5.1| , we discuss the production cross section for polarized J/ ip at the Tevatron. In 
order to phenomenologically implement the production cross section for B mesons at the 
Tevatron, we introduce a simple two-parameter flt procedure. In section P, we describe 
the relevant NRQCD matrix elements which we use for this analysis and also discuss the 
various sources of input errors for our estimates. Following this, in section |^, we give our 
detailed numerical estimates for the polarized cross section and predictions for the helicity 
parameter, a. We also discuss the relevance/influence of the various theoretical input 
errors to our predictions. The differential cross section for J/ip and ip' production from B 
decays and the corresponding polarization parameter a are displayed and compared with 
current experimental data |Q. Since, the data includes feed-down channels from excited 
charmonium states, we account for this in our analysis to derive the polarization cross 
sections. Finally in section ||, we conclude with a summary of the results and comment on 
further possible improvements to the precision of our calculation. In the appendix of this 
paper, we have tabulated all the relevant matrix elements and their sources and give the 
values which we use in our analysis. 



2 Basic formalism 



In |10|, it was shown that effects of the lower momentum scales of the order, rriQ -v, mq -v"^ 
and Aqcd can be factored into matrix elements that are accessible only via non-perturbative 
techniques or from experiments. On the other hand, the short distance contributions that 
occur on scales larger than the heavy quark mass, can be calculated within perturbative 
QCD. A matching prescription is required to identify the perturbatively calculated short- 
distance part and the non-perturbative NRQCD matrix elements. In the following, we 
recollect this matching procedure for polarized quarkonium production which we shall 
later use for our calculation. This is followed by a basic description of the matrix elements 
and their scaling properties. 
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2.1 The matching procedure 

Let us consider an inclusive production of a quarkonium state H with momentum k and 
helicity A via a parton level decay process of the type / H{k, A) + X. The semi-inclusive 
decay width is given as 



J2 dr{I - Hik, A) + X) = ^j^J^ 5^(27r)^5W(p - k - Px) |r,^H(.,A)+xP , (1) 



where E and p are energy and momentum of the decaying particle, Eh is the energy of the 
quarkonium and the sum over X includes the phase space integration for the additionally 
produced particles. On the other hand, from the NRQCD factorization theorem, the decay 
width in (|I|) can be factorized into short- distance coefficients and long-distance matrix 
elements of local four-quark operators. Formally, 

1 d^k 

J2 dr{I - H{k, A) + X) = — 3^^ C--(P^ ^) X (^^n^) ■ (2) 

In (^, the short- distance coefficients Cmn and n denote some quantum numbers of the 
various states) depend only on kinematical quantities such as momenta and masses of the 
involved particles. They include effects of distances of the order l/mg and smaller, where 
tuq is the mass of the quarks from which the quarkonium H is built. The matrix elements, 

H(X] 

{Omn ), are expectation values of local four-quark operators, sandwiched between vacuum 
states, (0| . . . |0) . These cannot be calculated perturbatively, but are extracted from 
experiment or from lattice calculations. 

As a passing remark, we note that if the decaying particle is a hadron, then the parton 
level decay width in (0) must be folded with a suitable distribution function for the parton 
in the initial hadronic state. In this case, the factorization approximation requires the final 
quarkonium state to carry a large relative transverse momentum compared to Aqcd pl| , pi? |. 



The typical four-quark operators which are related to the long-distance matrix elements 
have the general structure 

= ^^IC'lxVHwx^lCn^ , (3) 

where ip and x are the heavy quark and antiquark non-relativistic field operators, re- 
spectively, and ICn and /C'|„ are products of spin and color matrices as well as covariant 
derivatives. Here, Vh{x) is a projection operator that projects onto the subspace of states 
that contains the quarkonium state H{X) and in addition soft hadronic final states denoted 
by Sx- These soft states are supposed to be light, due to the NRQCD cut-off requirement, 
i.e. their total energy has to be less than the NRQCD ultraviolet cut-off A to avoid double 
counting. Hence including them, 'Ph{\) is written as 



Vh(x) = 2^ \H{k = 0, A) + Sx) {H{k = 0,X) + Sx\. (4) 

Sx 
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The matching procedure between the complete theory and the NRQCD expression requires 
that the normahzation of the mesonic states in both the frameworks, i.e., in ([^) and (H) be 
the same. In this analysis, we follow the relativistic normalization procedure as suggested 
in The matching condition is given as 

J2i^1T)'6^'\p - k - px)V^,-,+x'Tl^c,^X = J2 ^rnn{V. k) {^J^ IC'IxVh WX^ ICn^P) • (5) 

X m,n 

To carry out the matching procedure explicitly, both the l.h.s. and the r.h.s. of (^) have 
to be expanded as a Taylor series in q and q'. The short-distance coefficients can then be 
simply identified by an order by order comparison in the coupling constant along with q 
and q'. 



2.2 Expansion and simplification of matrix elements 

As mentioned above, an expansion of the matrix elements is necessary for matching with 
the complete theory. In addition, by applying the symmetries of NRQCD the matrix 
elements are simplified and expressed in terms of standard matrix elements. The relative 
magnitudes of each of the matrix elements can be estimated, using velocity-scaling rules 
which we list here. 

In the case of J/ip, the independent matrix elements can be determined by simple tensor 
analysis, since it is a vector meson with J = 1. Therefore, the helicity label. A, transforms 
like a vector index in a spherical basis. This corresponds to choosing circular polarization 
vectors as the basis vectors for the polarization states of J/ip. A unitary transformation, 
given by the matrix 

/ -I/V2 -t/V2 \ 
4=[ 1, (6) 

V I/V2 -i/V2 / 

connects the two basis. Here, i runs from 1 to 3, whereas A takes the values +1, and — 1. 
In what follows, we list the matrix elements that appear in the calculation of ip production 
in b quark decays. Rotational symmetry as well as heavy quark spin symmetry allow them 
to be expressed in terms of standard matrix elements with well defined spectral states, 
^■^"•"^Lj. Details on the expansion and reduction of NRQCD matrix elements can be found 
in For the simplest matrix elements without any vector indices one finds 

{i^^xrj/^ix)X^i^) = Im^ioi/'^CSo)) , (7) 

{^^T^xVj/^ix)X^T'^^) = Im^iOi^^CSo)) , (8) 
for the color singlet and octet case, respectively. The factor Arric originates from the differ- 



ent normalization of states in []TU| and The remaining dimension 6 matrix elements 

can be reduced to 

(^t^T«X^j/^(A)xV^T» = le^U^m^iOi^^CS,)) , (10) 
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up to corrections of order f^. In case of the matrix elements with four vector indices we 
have 

(^t (-f^') a'T'^xPc-cc-c'X^ {'^^1 ^''^"'^) = 4ef e^"5'="me(08^/'^(=^Po)) , (12) 

which also receive corrections at order f ^. Other useful relations can be established on the 
basis of heavy quark spin symmetry. It relates the matrix elements with the same orbital 
angular momentum L and different total angular momentum J to each other. For instance, 
the P wave matrix elements at order f ^ are equal up to a multiplicity factor 

{O'J^CPj)) ^ {2J+l){Oi/^('Po)) . (13) 

Furthermore, the color singlet matrix elements can be related to the non-relativistic quarko- 
nium wave function, whose radial part is denoted by R^, evaluated at the origin. This can 
be achieved by means of the vacuum saturation approximation (VSA). 

(^V^X^j/^(A)xV^'V^) = e^ef^M^\R^\^ , (14) 

where due to the VSA an error of order O (ferric |-R^p) is induced. 

As mentioned earlier, the non-perturbative matrix elements relative importance is de- 
termined according to the velocity scaling rules. These velocity-scaling properties for the 
J/ip production matrix elements are summarized in table ^ 



matrix element 


f-scaling 


{Oi'^'Pj)) 

(oi^^epj)) 


^3 

^11 



Table 1: Summary of velocity scaling rules for NRQCD matrix elements 



3 The process h Jj^ + X 

Having given the basic formalism, we now proceed to describe J/'?/' production from h 
decays. The decay of a 6 quark is a weakly induced process and is described by the 
exchange of a IV boson, transforming the h into a c quark. The W subsequently decays 
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into a cqj pair, where the qj is a hght state and can be either a s or a (i quark. At the scale 
/i ~ mf,, after integrating out the W boson, the effective QCD corrected Hamiltonian 



induces the b ^ cc + qf transition. The relevant operators are 

Oi = c-fi,LcqfYLb , 

Os = c^,LT'^cqfYLT% , 



(15) 



(16) 
(17) 



with L = 1 — 75. Here, C[i](/i) and C[8](/x) are the effective Wilson couplings for the 
color singlet and the color octet operator, respectively. The Wilson coefficients should 
not be related with the short-distance coefficients of NRQCD, denoted by Cmn- To get a 
qualitative feeling for the relative strengths and their dependence on the factorization scale 
/i, we have shown the couplings in figure |I] for /i ~ . . .rrif,. Here, as{Mz) = 0.119 has 
been taken, which corresponds to ^qcif'^ = 93.14 MeV at LO. In particular, it has been 



o 




Figure 1: Dependence of the Wilson coefficients C[„](/i) on the scale /i. 
noted that the color singlet coefficient Chi exhibits a strong dependence on /i and even 
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vanishes near /i ~ rric at LO |^0[. This behavior which hints towards large higher order 
corrections can not be cured at NLO We will later allude to this problem, when we 
discuss the various uncertainties pertaining to our results. 



3.1 The decay width 

Applying the effective Hamiltonian (|I^) to the 6 — >■ cc + g/ decay, we calculate the matrix 
element, T at LO to be 

G 

+ qs] ■ c{ph,LT^c{p)qf{pfh^LT%{p,) } . (18) 

The four- momenta of the outgoing c and c quarks can be expressed as, p = \k + Lq and 
p = — Lq, where k is the total four-momentum of the cc-system and q is the relative 
three-momentum in the cc-rest frame. Lf is the Lorentz boost matrix that connects the 
two frames. Expanding T to linear order in q, we get 

+Cn ■ qf{pf)YLT%{pk) [2Lf {m,^^a'T-r] + leU^q^a'T'^r,) - k^^^T^r,] } . 

(19) 

The boost matrices, , also have to be expanded in q and to linear order are given to be 



ii = '5.. + kf-i)*'*i. (21) 



where the hats denote unit vectors. At LO, |Tp factorizes into a product of two rank two 
tensors: 

\%^c-e+,,\' = \Gl\V:,V,f\'W^^ (T« + TlS) ■ (22) 
W^y describes the transition of a 6 quark to g/ and has the simple form 

W^u = \- {{^f + mf)^^L{j,^ + mb)7,L) = 4 ■ p^p^ {S^^^u - ^^a^.Xu) , (23) 

with 

SfMauX = dfiaOuX — Qi^iuQaX + Q^iXQua , (24) 

and e^o-i/A being the usual antisymmetric Levi-Civita symbol. T^u'' refers to J/ip production 
in either color singlet (n = 1) or color octet (n = 8) channel and is obtained as 

• [2L'; (m.r/'tor^e' + ^6,n,r/'tg'V^e') - k^^^'] , (25) 
Tg) = l^abGfsy [2Lf [m.^^a'T'^r] + teiki^^q'a'T'^r]) - k'^^^T'^r]] 

■ [2^ {m.i^a^T'i' + ie,nj>v'h"'^''T'i') - k'^v'^T'C] ■ (26) 
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In the context of the PM, the tensor W^^, will be replaced by a more general hadronic 
tensor structure which includes a distribution function for the heavy b quark. This will be 
discussed in section ^. In the following, for clarity, we describe in some detail, our calcu- 
lation for obtaining the polarized decay spectrum. Contracting the two tensor structures 
T/i"'* and Wfj_^, one can identify six different non-relativistic four-quark operators which 
then will be matched to the corresponding NRQCD operators that have been presented in 
section ^ The contraction can be divided into two steps, i.e., (i) the contraction of the 
Minkowski indices and (ii) the contraction of the three-vector indices. In the first step 
only four quantities have to be calculated, since T^"'' consists of the structures I/^I/^ and 
kfj,ku (mixed structures like L^^ki, vanish by symmetry arguments) and W^i, consists of a 
symmetric part, proportional to S^^xu, and an antisymmetric part, proportional to e^^xu- 
The four quantities will be denoted by P with a superscript (s) for symmetric and (a) for 
antisymmetric, specifying which part of W^i, they come from. The symmetric terms are 



P. 



I] 4 ■ p"fplScr^j,XuLiL^j 

P(^) = 4 ■ p'^fplS^.x.k^k'' = 8{pfk){p,k) - A{p,pf)k^ . 



Similarly one gets for the antisymmetric part 



(a) 



4^ ■ p^p^t^^xu^L] = Aie^^xuk'^pmi] = Mp^Vk^ (-e^jmLxm) 



(27) 
(28) 

(29) 



pW = -4^.p^p^e.^,,Pr = 0. (30) 

Following this, we choose a reference frame in which the decaying b quark moves with 
arbitrary three-momentum, |p^|, and corresponding energy, Ei, = \/ml + Ip^P. The J/ip 

three-momentum is denoted by |fc| with energy, = ■sJm^ + \k\'^. The three- vectors, p^, 

and k, enclose an angle "i?. For this choice of reference frame, the above projectors (^) - 
(B0|) are evaluated to be 



P, 



Eb\k\ 
2mr 



+ 1- 



2m, 



Ip^I cos-i? kikj 



Eb\k\ 
2mr 



+ 1- 



2m, 



Ip^I cos-i? 



kiPbj + Pbikj + PbiPbj 



+ 4 [ml - Ei,E^ + \Pf,\\k\ COST?] 6, 



I] 5 



P('^ = 8{EbE^ - \pf,\\k\ cos -A {ml + EbE^ - jpj |fc| cost?) (2me)^ 



Pjj Aimc^ijfji 



Eb\k\ 
2mr 



1 - 



E^ 
2m, 



pj COST? 

"^m Pbn 



(31) 
(32) 
(33) 



Collecting together (^ - 12^) and (^ - 133|) the color singlet contribution is given to be 
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and for the octet we obtain 



In order to perform the matching procedure described earher, we need to insert the above 
results into ( P^ to get the squared matrix element, \Tb^cc+qf\^- Following the results of 
section |], we identify the short- distance coefficients Cmn by making use of the matching 
condition As stated before, the integration over the phase space of the additionally 
produced particles (which in our case is the g/ quark) has to be included into the sum over 
the hadronic rest X on the l.h.s. of (|^). Using the standard identity 



d^Pf 



d^Pf 



and performing the four-dimensional phase space integral over d'^pf, due to the presence 
of the S^^\pb — k — Pf) function, pf gets replaced by ph — k. Thus, the matching condition 
is obtained to be 

^2'K6[{ph-kf-m'j]e[{pb-k)o]%*_^^^,_^q^rb^cs+qf = ^Cmn{Pb,k){ilj^lC'l,x'Pj/i>{x)X^J^n^) ■ 



If 



(37) 

We begin with the color octet contributions whose short-distance coefficients we identify 
for the spectral states '^^'^^Lj of the cc pair. In the following, the three-vector indices on 
the l.h.s. have been suppressed. 



^Gl\V:,V^f\'Cl^6{{p,-kf 



m 



X 



1 2 

-EbE^ + \pi\\k\ cos?9) + {EbE^ - cos^9) 



(3J 



m. 



X i2 



Ri\Ph\, \k\,cos{})ki -phi Ri\Pb\, \k\,cos{})kj - pbj 



+ (ml — EfjE^ + \Ph\\k\ cos??) 6ij — 2mci R{\pf^\, \k\, cos ^)k 



Pbn 



(39) 



CspPo] = fGl\V:,V^f\'q,^S{{p, - kf - 



m. 



X 



R{\Pb\, \k\,cosi^)ki-pf,i R{\Pb\, \k\,cos^)kj - Pbj 
+ [ml - EbE^ + \Pb\\k\ COST?) {6knSip - SkpSin) 



jnp 



— 2mri 



R{\Pb\^ \k\,COS'd)km -Pbm {Smieknp - Smk(^l 



LnpJ 



m 



2 ' 



(40) 
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where the kinematic function 



Ri\Pb\, \k\,cos^) 



Eb\k\ 
2mr 



2mr 



Ip^I cost!} . 



(41) 



The color singlet short- distance coefficients Ci[^^'^^Lj] is most easily obtained from the 
corresponding octet coefficient by replacing the color matrices by unit matrices as well 
as changing the Wilson coefficient from C^s] to along with an overall factor of |. This 
also serves as a useful book-keeping device for our calculations. 

/^From table |l|, it is seen that the color singlet matrix elements with angular quantum 
numbers ^5*0 and ^Pq scale with relative to the baseline matrix element. Additionally, 
the color singlet production is suppressed relative to the color octet production. This 
follows from the comparison of their Wilson coefficients, whose squared ratio turns out to 
be C|]/Cp^j ~ 25, which can be estimated from figure |I]. Therefore, the contributions of 
these color singlet matrix elements are highly suppressed and one only needs to take into 
account the ^Si contribution for the singlet case. On the other hand, all the three octet 
matrix elements should be included, because they scale as f ^ relative to the dominant 
color singlet ^5*1 matrix element and are enhanced due to the larger Wilson coefficient C^s]- 
Thus, we have the singlet contribution 



X 



\k\,cosd)ki -pbi R{\Pb\, \k\, cos ^)kj - pbj 



+ [ml — EbE^ + cos§) 6ij — 2nici R{\pi\, cos§)k 



Pbn 



(42) 



In order to calculate the decay rate, we choose a reference frame where the J/ ip moves along 
the positive z axis. The h quark momentum vector is then most conveniently parameterized 
in spherical polar coordinates, where d is the polar and ip is the azimuthal angle. We have 
the unit vectors 



ki 



5i3 



Pbi = COS ip sin {}6ii + sin ip sin {}6i2 + cos "OS, 



(43) 
(44) 



Multiplying the short-distance coefficients C„[^'^"^^Lj] from (|3^- ^) along with the appro- 
priate matrix elements, we get the differential decay rate as in (Q). After contracting the 
short- with the long-distance part we are left with a triple differential decay width for a h 
quark, moving at arbitrary momentum \p^\ that decays into a, J/ip with helicity A. The 
individual matrix element contributions in this case are 



{<^Ao {Eb\k\ - E^\p^\ cos^f - \p^\^(\ - cosH) + X [-Eb\k\ + E^\p^\ cos?9] 

+ l[ml+ \p,\\l - cos^^) - EbE^ + \p,\\k\ cos^] } x {Ol^'^CS,)) , (45) 
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-EbE^ + \p,\\k\ cos^ -ml + ^^ {E,E^ - \p,\\k\ cost?)'] } x (Oi^^CSo)) , (46) 



E^Ef, 



GUv:,V,f\'CLmMPt - ky 



m. 



{<^Ao [2^ {Eb\k\ - E^\p,\ cost?)' - \pf{l - cos' t9)] + A [-Eb\k\ + E^\p,\ cost?] 

+ l[ml+ \pf{l - cos't?) - E,E^ + \p,\\k\ cost?] } x {Oi^'^CS,)) , (47) 



E^Eb 



GUv:,v,f\'cLmMpb - ky 



m, 



[5x0 [\\Pb\\^ - cos' t?) - ^ {E^\k\ - E^\p^\ cost?)'] + A \-E^\k\ + E^\p^\ cost?] 



+ 1 



[|Pfe|'(l - cos' t?) + 2m,' - 2EbE^ + 2|pfe||fc| cost? + (Eb|fc| - cost?)' 

X ^(Os'^^^'Po)) . (48) 



In the above, 1 denotes the corresponding terms that contribute to the spectrum but have 
no exphcit A dependence. As a consistency check for our results derived so far, we choose 
Ip^I = 0, which corresponds to the h quark rest frame. Integrating over cost? and 
we obtain for the color singlet and octet contributions to the polarized decay width 



ri(6^ J/t/>(A)+X) 
and 

13(6^ J/t/>(A)+X) 



967r 



^F\^c\ycf \ ^ -Or- 



mlrric 



X 



{{ml - 4m,')5,o + 4m^(l - A)) {Ol^^CS.y (49) 



{{ml - Aml)5xo + 4m^(l - A)) (^^"/'^(^Si)) + ml{Oi''f' ^ S^)) 

+ 3 ((m^ - Aml)5xo + m,' + 4m,'(l - A)) 4(08'/'^(^Po)) 
respectively. This agrees with the result of Fleming et al. |3l[] . 



(50) 



3.2 Excited quarkonium states 

In the previous section, we applied the NRQCD factorization formalism to J /ip production 
in h quark decays, but with relatively small modifications we can equally well apply our 
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calculation to other quarkonium states with J = 1; as for example ip' and Xci, which are 
25" and IP states of charmonium, respectively. 

For production from h quark decays, only the J /if) matrix elements have to be 
replaced by ip' matrix elements, whereas the short-distance coefficients are not affected by 
this; 

{Oi'^^'^'Lj)) {OfC'^'Lj)) . (51) 

Besides, the inclusive decay width to Xci requires minor modifications, and the formalism is 
similar to the calculation developed for J/ip. Xci, being a ^Pi state of charmonium, to lowest 
order receives contributions from the (Oi^^l^Pi)) and (Oq^^I^Si)) matrix elements, other 
contributions are down by at least v^. The fact that to lowest order, a color octet matrix 
element significantly contributes to the spectrum also explains the difficulties to describe 
P wave quarkonium production in the framework of quark potential models. However, in 
the framework of NRQCD these problems are not completely resolved, even at NLO, since 
one faces the task of describing the production of all three Xcj states with one set of matrix 
elements. This is so, because the matrix elements are related to each other by heavy quark 
spin symmetry. Usually, they are expressed in terms of Xco matrix elements: 



(o^^e5i))=(2j+i)(ore^i)) . 
(o^^ep,))=(2J+l)(or('^o)) . 

Explicitly, the contributions to the differential decay width for b — > xa + X are 



(52) 
(53) 



d'T,['Po] _ 1 



{ho [11^/(1 - cos^ ^)-2k ^^^1^1 " ^xlPbl cos^)'] + A [-E,\k\ + Ejp,\ cosi^] 
+ l[\p,\\l-cosH) + 2ml-2E,E^ + 2\p,\\k\cosi^ + ^^{Eh\k\~Ejp,\cosi^)'''j } 

X -L(Of°(3Po)) , (54) 



E^Eb- 



dk'^ 



12n^ 



{Eb\k\ - E^Ip.I cos^9)' - \PbWl - cos' ^9)J + A [-Eb\k\ + E^\p^\ cost?] 
+ l[ml + |p6p(l - cos'i?) - EbE^ + \pb\\k\ cos^9] | x {Of\^Si)) . (55) 



3.3 Feed-down channels: H 



Apart from direct J/ if) production which accounts for roughly 70% of the J/ip from B 
decays, there are contributions from feed-down channels. Following a few simplifying 
assumptions as described in |jl8|,Q, it is possible to incorporate the J/ip production from 
these feed-down channels. 
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The first step towards this is to calculate the momentum spectra for the excited char- 
monium states. In the case of ip' , because it is an S state the procedure is the same as 
for J/tp, and for the Xa production rate, the necessary modifications are very modest as 



stated in section |3]^. Next, we need to evaluate the production rate of J /ip from ip' and 
Xci decays. It is assumed that in the excited charmonium decays, the three-momentum is 
transfered completely to the J/^, i.e., the ip' and Xci differential production cross sections 
are simply multiplied by their experimental branching fraction to J/ip final states. The 
different helicity states are taken care of by additionally weighting the helicity dependent 
production rates for ifj'y^ and Xci{x) with probabilities P{H\ J/'ipx') where H = ip',Xci, 
which describe the transition of a ijj' or Xci in helicity state A to a J/^ in helicity state A', 
respectively. 

ip' dominantly decays hadronically into J/ip, and since no spin flips are observed the 
polarization is unchanged by this process [|l8l. Thus we have -P(V'a ~^ J/i^x') = ^xx'- 
For the Xci state, the situation is somewhat different, because it decays radiatively into 
J/ip. The transition probabilities have been determined to be [^: P{Xci(o) ~^ J/i^o) = 0, 
^(Xci{±i) J/ipo) = |, ^'(Xci(o) J/^±i) = 1 and P(Xci(o) J/i^±i) = i Generically, 
the J/ip production through feed-down channels can be summarized as follows 

dT{B ^Hx^ J/tpy) = dT{B ^ H^) ■ Br{H ^ J/^) ■ P{Hx ^ J/tpy) , (56) 

where the inclusive parts X have not been noted in the transitions. Apart from the above 
mentioned feed-down channels, there can also be radiative transitions from Xco and Xc2 
which have been observed. Their production rates or branching fractions to J/ip are small 
compared to the ones of ip' and xa and hence have been neglected in our analysis. 



4 The hadronic decay B ^ J/ip -\- X 

The results derived in section ^ describe J/ip, or more generally, charmonium production 
from a free b quark decay. In the b rest frame, J/ip is produced with fixed momentum, 
because the kinematic implications of soft gluon emission in the J/ip production process are 
neglected. Therefore, at leading order in a^, the momentum distribution of the charmonium 
state results only from the Fermi motion of the b quark in the B meson. To incorporate 



these bound state effects of the B meson we adopt the PM as introduced in p9|,^p4[ 
An application to semi-leptonic B decays was first proposed in [^|. 



In all our calculations so far, the b quark occurs exclusively in the b —> qj transition. 



described by the tensor W^i, in (|23D . Introducing light-cone dominance as in p3i^, it 
is possible to relate the hadronic B — * Xf transition to the heavy quark parton distribu- 
tion function (PDF) f{x). In contrast to semi-leptonic decays, in inclusive charmonium 
production the momentum transfer is fixed, if one neglects the kinematics of soft gluon 
emission of the final charmonium state. Due to the on-shell condition for the charmonium, 
we have = k"^ = ^ 10 GeV^ which justifies the light-cone dominance assumption. 

At the computational level, the W^^ tensor is modified in two ways; (i) the quark 
momentum pb is replaced by the fraction of the B meson momentum xps, and (ii) the 
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entire partonic structure is folded with the PDF f{x) for the heavy b quark. Thus, we 
obtain 



W, 



flU 



4:{Sf,a,yX - i(^t^au\) / dxf{x)p\{xpB - kye[{xpB - k)Q]6[{xpB - kf - m 
Jo 



with the sign function 



eix 



+1, 
-1, 



x>0 
X <0 



(57) 



(58) 



This modification leaves the T^u^ tensors unchanged, they remain as in (^) and (|26|) . The 
distribution function dependence on the single scaling variable, x, is a consequence of light- 
cone dominance. In this framework, the distribution function is obtained as the Fourier 
transform of the reduced bi-local matrix element at light-like separations, hence 



fix) 



(59) 



as shown in [^] . 

In an infinite momentum frame, the distribution function is exactly the fragmentation 
function for a high energy b quark to fragment into a B meson Hence, the Peterson 



functional form can be adopted as a distribution function for the heavy b quark inside 
the 5 meson, with 

Here, f{x), is a one parameter function with a free parameter, ep, while A^^ is a normal- 
ization factor defined such that, 



dxf{x) 



(61) 



i.e., with unit probability there is a 6 quark in the B meson. In the parameter range that 
is usually chosen, ep ~ 10^^ . . . 10~^, f{x) peaks at large values of x; a behavior that has 
been determined elsewhere p7| , ^ . For completeness, in figure 0, the x dependence of f{x) 
is shown for four different values of ep. 

As a useful consequence of the PM, the b quark mass, rrib, is replaced by the B meson 
mass Mb removing the uncertainty in the quark mass. On the other hand, the parameter 
Ep now carries quite a large uncertainty, and in some form, we have traded one uncertainty 
for another. For comparison with other calculations, it can be useful to define an effective 

nib = {x)Mb , 



b quark mass |^ | 



where 



dx xf{x) 



(62) 



(63) 
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Figure 2: Peterson function for ep = 0.004, 0.008, 0.012 and 0.018. 



To illustrate the dependence of (x) and on ep, we give their values for different choices 
of ep in table | setting Mp = 5.279 GeV. 



As expected, the results of - ^8]) get modified due to the new hadronic tensor W^,y 
in (^). The integration over the scaling variable x is straight forward and therefore can 
be performed immediately. The generic integral / under consideration is of the form 



dx G{x)9[{xpB — k)Q]6[{xpB — kY — 



(64) 



where G{x) is an arbitrary function of the integration variable x. Note the presence of 
the 9 function which ensures that the final state qf quark has positive energy^]. This is 
required, because in the inclusive approach qf has to hadronize, finally giving a hadronic 
state Xf which of course has to have positive energy. 

To solve the integral, we introduce a new integration variable according to 



[xps — kY — m'j 



Ml 



{pBkf 



+ 



Ml-m) ^ 
Ml 



(65) 



^The function thus effectively replaces the e function in (57). 
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ep 


(x) 


rrih in GeV 


0.004 


0.85 


4.48 


0.008 


0.81 


4.29 


0.012 


0.78 


4.16 


0.018 


0.76 


4.02 



Table 2: Effective b quark; mass and expectation value (x) of f{x) for various ep and 
Mb = 5.279 GeV. 



with 

This translates (0) into an integral over z and 

Ml) 1 
/=/ dzG{z + zo)9[{z + zo)pb{o) - ko]-j-p2 — iS[z - zi] + 6[z + zi]) , (67) 

^2(0) ZIVl^Zi 

with 

zo = Psk/Ml , (68) 



^1 = ^{PBky - Ml{Ml - m))/Mj, . (69) 

Here, we have used the standard identities for the 5 distribution. The second 5 function in 
(|67|) does not contribute due to the 9 function which restricts contributions of the integral 
to the argument of the first 5 function. Therefore, ( p^ ) can finally be rewritten as 

1 

2M|2i 
where we have defined 

a;+ = zo + ^i = ^ [psk + ,J {pskf - Ml{Ml - m))) . (71) 

The two 6 functions express the fact that the scaling variable x^, since it can be interpreted 
as the b quark's momentum fraction within the B meson, is only allowed to vary between 
and 1. 

Inserting the above expressions into (^ - ^) and together with the tensor (|57|) , we 
obtain the following expressions for the decay width, sorted by their production mechanism. 

[5xox+ 2^(^B|fc| - E^\pb\ COS i3f - |pbP(1 - cos^^?) + \[-EB\k\ + E^\pb\cosi9] 
+ 1 [x+Ml + x+\pb\\1 - cosH) - EBE^ + \pj,\\k\ cos^] }x {Ol^^CS,)) . (72) 



18 



1 



Tfl 1 



{l \-EBE^ + \pB\\k\cosd-xMl + ^x+ {EBE^-\pB\\k\cos^f'\ } x (08^/'^(i5o)) 



(73) 



E,i,E 



^ r"^ IT/*T/ |2/^2 "^c 1 
?B Zi 



{'5Aoa;+ 2^(^B|fc| - E^\pb\ COS i3f - - cos^^?) + A[-^B|fc| + E^lp^l cosi?] 

+ 1 [x+M| + - cos' 7?) - EbE^ + \pj,\\k\ cos^] } x {Oi^^CS,)) . (74) 



dk^ 



247r2 



'AO 



- cos'^) - 2^x+ (Eij|fc| - E^\pb\ cosi^f + X [-EB\k\ + E^\pj^\ cos^9] 



+ 1 



x+|pb|'(1 - cos'^9) + 2x+M| - 2EbE^ + 2|ps| |fc| cos^9 

+ 2^^+ (^i^l^l - EM cos^)'] } X -i,(03^/^(3po)) . (75) 



The distribution function f{x), being defined in an infinite momentum frame, as well as 
the incoherence assumption restrict the PM application to B mesons with large, if not 
infinite, three-momenta. To circumvent this problem, in the case of unpolarized decay, a 
Lorentz invariant quantity, E ■ dT, is usually constructed and subsequently evaluated in an 
arbitrary reference frame. This strategy was adopted in the calculation of semi-leptonic B 
13) and also in where the PM was first applied to inclusive hadronic 



decays (e.g. HS 



B decays. This enabled the authors to employ their calculations to B decays at the CLEO 
experiment where B mesons are produced almost at rest. However, this is not possible in 
the case of polarized production cross sections, since they are frame dependent. Therefore 
one has to go to a large momentum frame to fulfill the requirements of the PM. This brings 
us to the description and apphcation of our results to the Tevatron. 



5 B ^ J/t/j + X at the Tevatron 

At the Tevatron, B mesons are not produced with fixed momentum as in T decays 
(e.g., at CLEO), but in fragmentation mode. Therefore, one has to deal with a mo- 
mentum distribution of the B mesons, expressed by the pt and Vb {B rapidity) dependent 
double differential B production cross section. The double differential B cross section 
d'^a{pp —>■ B + X)/{dpTdyB) has to be folded with the normalized semi-inclusive differ- 
ential decay spectrum, 1/Tb ■ d^r{B ^ J/ip + X')/dk^, to obtain the desired differential 
production cross section for J/ip. 
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The CDF collaboration at Tevatron has already measured J / ip and ip' production from 
B decays []50|,pl[] as well as their polarization [^. For the latter we are not aware of any 
theoretical predictions that take into account bound state effects. Our results presented 
here (with some kinematic modifications) are directly applicable to the Tevatron exper- 
imental setup. We proceed to discuss the required kinematics to suit the Tevatron and 
make a comparison of our predictions with the available data. 

So far the absolute values of the three-momenta \Pb\i 1^1 the polar angle d between 
Pb and k have formed the set of kinematic variables in our calculations. For these, we will 
trade with a different set: the transverse momenta pt and kx-, the rapidities ub and y and 
the azimuthal angle 0. px and ub belong to the B meson, whereas fcr, y and describe 
the kinematics of J /ip. Following this, the relations between the old and new variables are 
given to be 



|fc| =(4 + sinh2y(4 + Mj))'/' , 
\Pb\ = {pt + sinh^VBiPT + MDY'^ , 

cos •§ ^jj^jj — I ^sinh k"^ + M| ■ sinh yB\JpT + ~ kxPr cos 0^ 
Here, the rapidities are defined to be 

'E + kn 



(76) 



y 



( Eb+P\\ \ 
and = - In — , 

\Eb-P\\J 



(77) 



with k\\ and p\\ being the momentum components parallel to the beam line. The Jacobian 
determinant for the coordinate transformation of the ip variables is 



d{kx-, ky, kz 



kTJk^ + Mlcoshy . 



(78) 



5.1 Cross section and the a parameter 

In order to make a comparison with data, we have to take in to account the relevant 
experimental constraints and cuts. In the case of J/ ip and ip' production, the cross section 
is rapidity integrated over the region \y\ < 0.6 0,|51|,0. For direct J/ip production the 
cross section is expressed as 



da{J/ip) 
dkr 



dpT / dyB 



d^ajB) 
dprdyB 



dy 



d{kx, ky, kz 



dikT,y, 



TndP 



{pT,yB,kT,y,(i)) ■ (79) 



As already discussed in section |3.3| , apart from the direct production also the feed-down 
channels contribute to J/ip production. These can be incorporated under the assumptions 
made in section |3.cj| and hence the decay rates to J/ip just have to be summed for the 
direct and the feed-down channels to obtain the final result. Using (|79D, we can evaluate 
the usual polarization parameter, a, which is experimentally accessible with the help of 
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a fit to the angular distribution in the di-Iepton decays of J/ip. The angular differential 
decay spectrum has the following form 



dcos6 



(xl + a cos^ e , (80) 



where the angle 6 is defined in the ip rest frame in which the z axis is aligned with the 
direction of motion of the J/ip in the B rest frame. 

Theoretically, a is expressed as the ratio of linear combinations of the helicity produc- 
tion rates for J/ip- Given the decay of a longitudinally polarized vector particle (helicity 
A = 0) and for a transversely polarized state (helicity A = ±), using (^) one immediately 
obtains 

0-+ + - 2(To X 



a 



Cr+ + CT- + 2(70 

The helicity production cross sections a\ are the ones obtained from ([79|), which makes a 
a function of kx- 



5.2 Present status 

With respect to the Tevatron, there are two ways of employing the differential B production 
cross section in our calculation; either theoretically or experimentally. One approach, 
is to calculate the cross section by means of the QCD improved PM. In this case, the 
partonic b quark production cross sections have to be calculated in perturbative QCD and 
subsequently have to be folded with (i) the non-perturbative PDFs of the incoming hadrons 
and (ii) the fragmentation function, describing the hadronization process of the b quark. 
In the second approach, one can directly use the measured B production cross section at 
the Tevatron c.m.s. energy y/s = 1.8 TeV. 

It is known for quite a number of years that the theoretical description of the b pro- 
duction process fails to reproduce the experimental data (see e.g. |§2|). The shape of the 
spectrum comes out as desired, but the normalization usually falls short by a factor of 2 
or more (e.g. An agreement with the data can be achieved if the involved pa- 

rameters like the factorization scale, /i, Aqcd and mh are driven to rather extreme values. 
Recently, there has been a suggestion that the fragmentation function ansatz which has 



been used most frequently might not be appropriate |5^. However, in a next-to- leading 
order calculation of B meson production in pp collisions, it was shown that it is possible to 
accommodate the data within experimental error bars without fine tuning of the relevant 



parameters In this case, the b fragmentation function has been fitted to LEP data 
and the differences in the scales between LEP and CDF data is accounted for by the usual 
evolution equations. 

At this stage, it is not clear if the corrections due to the b fragmentation are responsi- 
ble for the disagreements between theoretical predictions and experimental data. It may 
also be that the theory for the b production mechanism is incomplete. There has been an 
attempt to account for the discrepancy involving physics beyond the standard model ^6 . 
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In the present analysis, we adopt the second strategy, i.e., to use an experimental 
fit. We are motivated to this choice, because, the main concern of this work is not the 
dynamical mechanism of B meson production in pp collisions, but their subsequent decay 
into quarkonium states. A strong argument for this being, the available data have been 
extracted using exclusive —> J/'il! + decays with the CDF detector at the Tevatron 
5^,0. The same group has also measured the J/ip and ip' production cross section |^| 
and their polarization Therefore systematic errors in the analysis should be of less 
significance and the comparison between theory and experiment will be more meaningful. 
The problems due to the lack of experimental data for the double differential B production 
cross section will be discussed later in this section. To this end, we first derive a useful 
algebraic fit to the production cross section in the following subsection. 



5.3 A two parameter fit 

It has been shown that the differential production cross section with respect to the B 
transverse momentum, px-, exhibits a simple power law behavior with a power in the 
range of —3 and —5. We therefore choose a power law ansatz of the type 



da 

dpT 



A-pV, 



f82) 



with A and n being the two fit parameters. Here, we minimize the relative rather than the 
absolute deviation square. This is done to make the functional form for the fit to correctly 
reproduce the cross section for both low and high pt- The reason being, at high p^, the 
absolute value of the cross section becomes rather small because of the power law behavior. 

To estimate the uncertainty due to the experimental B production cross section, we 
apply the fit procedure not only to the central values, but also to the cross section ± sta- 
tistical errors, i.e., we fit the "upper/lower ends of the error bars". In table |^, the results 
of the three fits including their values are given, and these are also displayed in figure 

I- 





central 


upper 


lower 


A 


2.4635 ■ 10^ 


2.8084-10*^ 


2.1363- 10^ 


n 


-4.0049 


-3.9984 


-4.0187 


Xdof 


0.166/2 


0.437/2 


0.125/2 



Table 3: Fit results of the model A ■ p^ io the differential B^ production cross section 
da/dpT for central, upper and lower experimental values. 

However, given the kinematics relevant to the Tevatron, we note that the fit for the 
experimental cross section as an input distribution is only available in rapidity integrated 
form. On the other hand, one needs to convolute the double differential B production cross 
section with the normalized differential decay spectrum for B J/ip + X with respect 
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Figure 3: Experimental data and fit results for differential production cross section 
da/dpT- 



to transversal momentum pt, as well as rapidity ys (see equation (^)). In the following, 
we assume the double differential cross section to be only weakly dependent on i/b in 
the dominantly contributing ys range. Later we see that this assumption gets justified 
for our setup. The error in the numerical result of the differential J/tp production cross 
section due to this treatment can be estimated by including a purely phenomenological 
ys dependence of the B production cross section, based on theoretical calculations. We 
have tried to reproduce the rapidity dependence of the double differential cross section 
d'^a{pp B -\- X) / {dpTdyB) given in ||5^ for small values of ys and have chosen the 
normalization such that the central experimental values of da/dpx in were recovered 
when integrating over ?/_b in the range —1 to +1. In a wide pt range (10 ... 20 GeV), we 
note that, d^a / {dpTdys) reduces by roughly 15% between \yB\ = and = 1 and so the 
factorization assumption made here seems to be reasonable. A particularly simple form 
for the double differential cross section can be taken to be of the form 



d'^a 
dprdy 



B 



{B - Cy^S) ■ P^T 



^3) 
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Note that the function in ( p3D has to be even in y^. To satisfy this, the simplest choice 



is m = 2, which, for not too large values of i/b, reproduces the shape of the curve in [|55 | 
sufficiently well. For this choice, we have a few constraints which are to be satisfied; B and 
C have to be chosen such that the 15% reduction between \yB\ = and 1 as mentioned 
here is accounted for. Furthermore, the parameter A has to be recovered when integrating 
over the yB interval [—1, +1]. These conditions together fix the parameters unambiguously 
to be 

where Xy ~ 15%. We use this result to estimate the error in our fits for the production 
cross section when we assume yB independence. 

5.4 Cumulative cross sections from all h Hadrons 

The B cross section that we applied in the last section refers only to B^ production, 
whereas & B ^ ip + X decay could involve any B meson type and even A;, baryons. 



Analysis of h quark fragmentation fractions are available from LEP data and from the 
CDF collaboration (Tevatron) [^,^ and we use them to include contributions from 5'', 
i?g and Af, baryons. As a passing remark, we note that the number of produced B^ mesons 



is too small to give a sizeable contribution . We make the simplifying assumption that 



the shape of the production cross section is the same for all h flavored hadrons and hence 
we multiply the B^ differential cross section by a factor Fb to include the contributions 
from the other h hadron types. This is certainly a good approximation for 5°, but ignores 
the mass differences for 5° (Am ~ 100 MeV) and A;, (Am ~ 340 MeV). Finally, the net 
production rate has to be multiplied by 2, because at the quark level both h as well as the 
charge conjugated decays can give rise to charmonium states. 

Apart from the fragmentation fraction Br{h B^), which at the Tevatron is 
measured to be 2a above the LEP result, the extracted values are compatible within one 
la. We choose to use the Tevatron data |^ as central values for our analysis, due to the 
possible differences of b fragmentation in e~^e~ and pp collisions. The experimental data 
as well as the combined scaling factor, Fb, that we finally use to incorporate the other b 
hadron contributions are summarized in table 0. 



6 NRQCD matrix elements 

The most crucial input parameters in the entire calculation are the NRQCD matrix ele- 
ments, because their values influence the theoretical predictions significantly. Numerical 
values of these matrix elements for J/ip, ip' and Xcj production that have been published 
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fragmentation fraction 


CDF notation 


CDF value 


Br{b B+) 


fu 


0.375 ±0.015 


Br{h B°) 


fd 


0.375 ±0.015 


Br{b 5°) 


fs 


0.160 ±0.025 


Br{b -> A,) 


/baryon 


0.090 ±0.028 




2.67±0.11 



Table 4: Experimental results on b fragmentation fractions at the Tevatron and the 
combined scaling factor Fb- 



in the literature are summarized in the appendix (see section |A.1|) . They have been de- 
termined mostly in unpolarized prompt production at hadron colliders, but also at ep 
colliders and in fixed target experiments. Apart from the ^Pj element, all the other matrix 
elements considered here have to be positive [^. In the following, we quickly recapitulate 
the extraction of the color singlet and octet matrix elements from various phenomenological 
data which we shall use. 



6.1 Color singlet matrix elements 

Leading color singlet matrix elements can be obtained in various ways. A popular approach 
is to calculate the quarkonium wave function within a quark potential model, adopting the 
QCD inspired Buchmiiller-Tye potential Alternatively, experimental decay rates of 
quarkonium states can be used to obtain numerical values for the color singlet matrix 
elements. For S states, the partial decay width, r{il) e^e~), is expressible in terms of 
the color singlet matrix elements up to higher corrections in the expansion [^]; and 
hence, r{ilj — ^ e~^e~) oc {Of{^Si)). Using experimental data from [Q, we find for J/ip the 
values 

{O(^^CSi))^0) = (0.77 ± 0.06) GeV^ and {0(^'^CSi))^i) = (1.37 ± 0.10) GeV^ , (86) 

where the results are displayed without (subscript (0)) and with QCD corrections (subscript 
(1)), respectively. The error corresponds only to uncertainties in data. Here, we have 
chosen as{2'mc) = 0.26 and Mj/^ = 3.097 GeV. Evaluating the same formula for the ip', 
we obtain 



(Cf (=^5i))(o) = (0.44 ±0.04) GeV^ and (Cf (^5i))(i) = (0.78 ± 0.07) GeV' 



(87) 



where M^/ = 3.686 GeV has been adopted. For M^i = 2mc = 3 GeV, these values reduce 
by roughly one third. 

For P wave charmonium states, the electromagnetic decay Xcj — >■ 77 for J = and 
J = 2 are suitable to fix the leading color singlet matrix elements ((9f''"^('^Pj)), which are 
related to each other through (^3]). For the Xco decay, the precision of the available data 
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is lower than for the Xc2 decay [Q. Therefore, the latter one is adopted more often to 
fix the corresponding (O^^^I^Pq)) matrix element. Including QCD corrections to lowest 
order, the decay width can then be expressed in terms of the leading color singlet matrix 
element fl^. In this case, we have, T{xc2 — > 77) oc {Oi^'^^Pq)). Numerically, we find 



(C»f°(3Po))(o) = (3.9 ±1.4) ■ 10"^ GeV^ and (Of °(3Po))(i) = (7.0 ±2.5) ■ 10"=^ GeV^ (88) 

for the tree level and the QCD improved calculation, respectively. Here M^^^ = 3.556 GeV 
was chosen, M^^^ = 2mc = 3 GeV results in a 50% reduction of these values. 

Alternatively, the decays of Xco and Xc2 into light hadrons have been used to fix the 
above matrix elements |Tl| with results that are compatible with the ones obtained above. 

A third possibility which also is adopted quite often, is to calculate the matrix elements 
on the lattice [^. In principle all the three methods give matrix element of comparable 
size, in particular, if one chooses the mean value of the leading order and the QCD improved 
result in the experimental analysis. 

6.2 Color octet elements 

In the case of J/tp and ip' prompt production data from fixed target as well as collider 
experiments have been used for extracting the color octet elements. Due to the different 
krp dependence of the matrix elements contributions, it is possible to fit their values to the 
unpolarized differential prompt ip production cross section, da/dhx- At large /c^, due to 
gluon fragmentation dominance, as described in the beginning, the ^5*1 matrix elements 
become the dominant source for ip production. On the other hand, the gluon fusion process 
gives rise to contributions that are proportional to the ^5*0 and the ^Pq color octet matrix 
elements. Since these two contributions exhibit a very similar dependence, the matrix 
elements cannot be fixed individually in the fit procedure. Instead, it was suggested to fix 
a linear combination of both p7|, |6^] , which is conventionally denoted by 

Mf = {OtCSo)) + —,{OtCPo)) , (89) 

where r is empirically determined and lies in the range 3 . . . 3.5 for hadropro duct ion. 

The fitting procedure for the NRQCD cross section for hadroproduction, as well as for 
electro- and photoproduction involves PDFs as theoretical input of the calculation. Hence, 
a theoretical uncertainty is introduced, because for different PDF sets the values of the 
matrix elements differ quite severely (see e.g. [|69| ). 

There have also been efforts to individually fix the ^5*0 and ^Pq color octet matrix 
elements for ip production by fitting their values to the ip momentum spectrum from B 
decays for the spectra as measured by the CLEO collaboration []7D[. However, at 



the high momentum end of the spectrum, as was noted in [^, one encounters difficulties 



because of effects of resonant two- and three-particle decay channels in the spectrum. In the 
context of the Tevatron, since only the absolute normalization of the spectrum is of primary 



importance, it is best to fit for the branching fraction. This procedure was first used in p6 
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and the matrix elements were fixed to the inclusive branching fractions, Br{B ip + X) 
for the CLEO data. In our analysis, we adopt this procedure as mentioned in the following 
subsection. 

6.3 Prescription for analysis 

All the difficulties mentioned above, make precise theoretical predictions involving the 
matrix elements to be very hard. A possibility to constrain especially, the {0^{}Sq)) 
and the (Og (^Pq)) niatrix elements more severely is the calculation of the polarization 
parameter a. The primary reason is a being a ratio of cross sections is less susceptible 
to theoretical uncertainties from other sources, and we shall discuss this in section ^ 
Concerning the matrix elements, we will incorporate the theoretical uncertainties, which 
have their origin in higher order corrections and the fit procedure. This we do by letting 
the matrix elements vary in ranges which cover almost all values given in the literature. 



These ranges, are summarized in appendix |A.1| at the end of each table. On the other 



hand, to avoid the predictions becoming too loose, we will impose three conditions on the 
values of the matrix elements. 

1. All observables are restricted to physical values, i.e., for instance the momentum 
spectra, dV /d\k\, have to be positive for all values of 

2. The numerical values of (08(^5*0)) and (Og (^Pq)) are required to give Mf that lies 
within the determined range. 

3. The resulting branching fractions, Br{B ip + X), have to match with experi- 
mental data |]70||5^] within error bars. Hence, the variation of the matrix elements 
within their (rather large) errors is not independent anymore. Furthermore, ep, is 
constrained by this restriction, since its variation influences the branching ratio. 

The above conditions do allow us to constrain the theoretical uncertainties to a reasonable 
measure. 

In order to have a quantitative feeling for the influence of the B cross section on the 
J/ip production cross section (flgure ^) and the contributions of the individual NRQCD 
matrix elements to the J/ip production cross section (flgure ^), we flx a set of input matrix 
elements. Here, for these two flgures, we have used the following numerical matrix element 
values as input for calculating the J/ip production cross section and will refer to them as 
the standard set: 

{Oi^^CSi)) = 1.1 GeV=^ , 

{Oi^^CSi)) = 0.0075 GeV^ , (90) 
{Oi^'^CSo)) = 0.055 GeV^' , 
{Oi^^CPo)) = -0.006 GeV^ , 
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along with Sp = 0.012. These values correspond to M'r^'^ = 0.046 GeV^. We remind 
that our choice in (pOD is different from that made in [^,|29|, since we are using a generic 
representative value taken from the allowed range as shown in appendix |A.1| . The values 
are required to reproduce the experimental branching ratio, Br{B ^ J/tp + X) = (0.80 ± 
0.08)% taken from [|70|. We reiterate and clarify that our terminology "standard" is to be 
understood only in the context of the values in (|90|) which are used to reproduce figures ^ 
and ^ 

6.4 Other input parameters 
6.4.1 Quark masses 

The b quark mass, which is not very accurately known, does not appear in our calculations. 
This is because, by virtue of having adopted the PM, the quark mass is replaced by the 
corresponding hadron masses which are known to high accuracy |^2|. Nevertheless, an 



effective b quark mass, {mi,), can be defined in the framework of the PM as has been done 
in section ^ Hence, a variation of ep, is in some sense equivalent to that of mb. 

The value of the c quark mass is strongly correlated with the values of the matrix 
elements. We have chosen rric = 1.5 GeV and do not vary it independently from the 
matrix elements, because a variation of rric by ±100 MeV is mostly taken into account 
in the uncertainties of the matrix elements. So, the influence of the uncertainties is not 
analyzed separately. This value is also used in most standard publications that have 
performed fits of the matrix elements. The alternative choice, = is of course 

also appropriate in the case of J/ip, since it is only slightly heavier than 3 GeV. But this 
choice is unsuitable for the excited states like ip' and Xcj- For the latter charmonium 
states, it would severely reduce the branching fraction Br{B —>■ + X), if one does not 
perform a new extraction of the matrix elements with rric ~ 1-8 GeV. Another argument 
for the use of nic = 1.5 GeV, is that it is a short- distance parameter in the calculation 
whereas bound state effects which are responsible for the mass differences of the various 
charmonium states are long-distance effects and thus have to be included in the matrix 
elements. 

The light quark mass mj is in general set to zero in the numerical calculations. To 
estimate its influence on the final results we let it vary in the range mj = . . . 150 MeV. 
This has no significant impact on our results. 

6.4.2 Peterson parameter 

Along with the matrix elements, the distribution function parameter ep, is varied such 
that the theoretical branching fraction Br{B ip + X) coincides with the measured value 
within la, hence, this couples ep to the values of the matrix elements. The standard value 
for ep has been determined in fragmentation processes of high energy b quarks (using LO 



QCD calculations) and ep = 0.006 ± 0.002 |]7T|. On the other hand, a more recent LO 



extraction has found a somewhat larger value of ep = 0.0126 [BH]. In the case of semi 
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leptonic B decays, relatively small values of ep are preferred [Q, whereas, an application 



to hadronic B decays show better results for larger values, ep = 0.008 . . . 0.012 

For our analysis, we let ep vary between 0.004 and 0.018, which covers the whole 
spectrum of its allowed numerical value. We note that the relation between the fragmen- 
tation and the distribution function is exact only in an infinite momentum frame for the 



B meson |^5[. Further, the wide range of values for ep also reflects the fact that we deal 
with a spectrum of B momenta which certainly is not inflnite. Additionally, the range 
of (preferred) values for ep could be due to the different momentum transfer involved in 
semi-leptonic and hadronic B decays. 

6.4.3 Wilson coefRcients 



As already known for many years now [^, the color singlet Wilson coefficient C[i] is 



not only small compared to the color octet coefficient C[s], but also strongly dependent 
on the factorization scale, /x. For a particular choice of fi in the range, < fi < rrih, 
the singlet coefficient even vanishes thereby, completely switching off the color singlet 
contribution to charmonium production in B decays. This problem is not fixed even at 
next-to-leading order, and more so, it even gives rise to a negative, i.e., unphysical decay 
rate iO,!!- Solut ions to this problem have been proposed, e.g., an alternative combined 



expansion in as and C[i]/C[8] which seems to stabilize the evolution, but this proposal 



lacks a solid theoretical basis [^. We estimate the uncertainties due to the factorization 
scale dependence of the Wilson coefficients by letting fi vary between its standard value of 
4.7 GeV |0,|2^,^ to 3 GeV, and study the influence on the differential cross section and 
the polarization parameter, a. 

6.4.4 Sundries 

Other numerical parameters that occur in our calculation are the Fermi coupling constant 



Gp, and the CKM matrix elements iVrfP. We have taken the standard values 16^ which 



rp, anu ine ^^ivivi mairix eiemenis \Vcf 
are 

Gp = 1.1664-10"^ GeV~^ . 

iKbl = 0.0402 ±0.0019 . (91) 

If the light quark mass mj is neglected, one can simply sum over the flavor qf = d, s and 
apply unitarity to obtain the approximate relation 

\Vcdf + \V,sf-l. (92) 

For the B mass we adopt Mp = 5.279 GeV, because the B^ and B^, B^ dominate in the 
h hadron admixture that is encountered at the Tevatron. We take the average lifetime to 
be, Tp = (1.564 ± 0.014) ■ 10^^^ s. Branching fractions, Br{B — + X), that are required 
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to constrain the set of matrix elements are taken from p2| : 

Br{B J/^(direct) + X) = (8.0 ± 0.8) ■ 10"^ . 

Br{B -^f + X) = (3.5 ± 0.5) ■ 10"^ . 
BriB Xci (direct) + X) = (3.7 ± 0.7) ■ 10"^ . (93) 

The branching ratios for the feed-down channels are taken to be 



Br{tp' J/ip + X) = (55 ± OJ70 
Brixci J/^ + 7) = (27.3 ± 1.6)% . (94) 

To normalize the cross sections according to data the charmonium branching fractions to 
are required and we take them to be 

BriJ/ij fi^f^-) = (5.88 ± 0.10)% . 
Br{ip' /i+/i") = (1-03 ± 0.35)% . (95) 

7 Numerical analysis and discussion 

In this section we will present the numerical results of our calculations and compare them 
with available experimental data; the unpolarized quarkonium production cross sections 
from B decays [^[5l[] and the a parameter [|^. This analysis also illustrates the depen- 



dence of our numerical results on the various input parameters and the associated errors 
that are involved. 

7.1 B production cross section 

The uncertainties of the B cross section, da/dpjn, have quite a large impact on the J/ip 
cross section and is of the order of 15 . . . 20%. In figure ^, we show these uncertainties for 
unpolarized J/ip production cross sections when da/dpT is varied between its lower and 
upper limit, as specified in table ^ The matrix elements and Peterson parameter were 
adjusted to the standard values which are given in section |0| . For J/ip production, one 
clearly observes that the predicted cross section agrees quite well with the CDF data |^ 



for the range, kx = 7 ... 12 GeV and is within the statistical errors. For the rest of the 
values, it clearly overestimates the data. Here, it is not a theoretical failure, but rather 
an artifact of the fit procedure for the B cross section. To evaluate our formula in the kx 
range 4 ... 25 GeV, we need to extrapolate the B cross section to pr values as low as 3 GeV 
and as high as 75 GeV. Presently, data exist only for px = 5 . . .20 GeV. At the lower end 
~ 5 GeV, mass effects of the b quark start to play an important role, causing the cross 
section to stop rising as steeply as in the higher p^ region. Also for very large pt, our fit 
seems to overestimate the data quite severely. In the intermediate interval, where only B 
mesons with transverse momentum from the measured range contribute, we have a good 
agreement. 
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Figure 4: Uncertainties in da/dkr for J/ip production due to the experimental B cross 
section. 

However, as we focus on a which is a ratio of cross sections, it is almost unchanged by 
the variation of the input momentum distribution for b hadrons and we can extrapolate 
to very small/large pt- This makes the prediction for a insensitive to the errors due to 
the fitting procedure. Numerically, a varies at the level of 10~^, hence the errors can 
be safely neglected. The uncertainties of parameters that influence the differential cross 
section normalization have also been included in figure ^ Among these are the scale factor 
Fb, whose error is specified in table ^, the branching fraction, Br{ip and the 

CKM matrix element Vcb which are given in section |6.4.4 

The numerical impact of neglecting the rapidity dependence of the B production cross 
section is noticeable at the level of less than 3% in the calculation of the J / ip cross section 
and less than 1% in the results for a. The reason for this rather weak influence is that the 
experimental cut on the J/ip rapidity {\y\ < 0.6) also imposes an upper limit on ys- A 
simple kinematic calculation shows that the difference between y and yB is limited by the 
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following expression 



\yB -y\< cosh . (96) 

Ikl + Ml^p'^ + Ml 



For COS0 = —1 and y = ?/max = 0.6, we find an absolute upper limit on i/b as a function 
of kT and pt- In our calculation this expression reaches its maximum for fc-r = 4 GeV and 
Pt ^ 6 GeV, giving a numerical value of \yB\ < 0.95 . As stated earlier, the dependence 
of the double differential B cross section on is rather weak in this region and hence 
does not lead to significant deviations if it is neglected. This justifies our procedure in all 
posterity. 



7.2 Dependence on matrix elements 

Both the differential cross section and a are strongly influenced by the variation of the 
matrix elements. For a, it is the main source of uncertainty, while all other errors, are 
canceled to a good accuracy, since a is a ratio of cross sections. We also remind that 
the variation of ep goes along with the variation of the matrix elements, because of the 



conditions that we imposed in section |0 



The direct J/ijj production cross section shows a variation of ±10% on the matrix 
elements. It should be noted that da/dkT is not sensitive to the individual values of the 
matrix elements, but only to the value of the resulting decay rate r{B J/tp + X) or, 
equivalently, to the branching fraction Br{B ^ J/ip + X). The maximum value of the 
cross section corresponds to a combination of matrix elements and ep, leading to the upper 
limit value of Br{B J/ip + X) = 0.88%, whereas the minimum value corresponds to the 
lower limit Br{B J/ip + X) = 0.72% (see sections and |6.4.4|) . 



In figure ^, we show the individual contributions of each matrix element to the unpolar- 
ized direct J/ip production cross section, da/dkT, for a standard set of matrix elements and 

as given in section |6]^. It has to be noted that for {O^^^ {^Pq)), we show the absolute 
value of the contribution, since it is negative (see section]^). As can be seen, the transverse 
momentum dependence to the different contributions is approximately the same. Thus 
this feature does not allow us to concentrate on any one matrix element contribution in a 
particular kinematic region. This has to be contrasted to prompt production, where the 
contributions of the matrix elements have different kx dependence. 

To quantitatively illustrate a's dependence on the various matrix elements, we give the 
helicity parameter for direct J / production, which we find to be 

'i-^ (<^(350) + 6^(03^/^(350) + fer''((^8'^''m)) + &;"''(08^^(^^o)) ■ 

To avoid the strong kx dependence of the cross section, we have normalized the color singlet 
coefficient in the denominator to unity. Note that the numerator in (pT]) is independent 
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Figure 5: Contributions of the different NRQCD matrix elements to the direct J/ip pro- 
duction cross section da/dkr- 

of {Og^^(^So)), because its contribution to the decay rate is hehcity independent and thus 
cancels out in the numerator. The coefficients ^-^ and b^^^^-' are given in table |^ as 
functions of fc-r for ep = 0.012 . Clearly, a is sensitive to the numerical values of the matrix 
elements and also depends on the short-distance coefficients. Apart from the normalization, 
due to identical short- distance coefficients, a variation of both the singlet and octet ^5*1 
matrix elements have a similar effect on a, i.e., with an increase of the matrix elements, a 
decreases (its absolute value increases). For {O'^^^ Sq)) , the situation is extremely simple, 
since it only contributes in the denominator, hence an increase of its numerical value leads 
to a larger value of a (reduction of the absolute value). A particularly strong dependence 
is exhibited in the case of {Og^^ {^Po)) , as it can even change the sign. In fact, for central 
values of the other matrix elements, it is possible to find a reasonable value of {Og^^ {^Pq)) 
which leads to a diverging polarization parameter. The reason being the short- distance 
coefficient in the denominator in (^) is so large that a small negative value of the matrix 
element can result in a vanishing denominator. Obviously such a value of {O'^^ {^Pq)) 



violates the constraints of section 6.2, because it corresponds to a negative value of the 
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kx in GeV 


ai 




ag 




"8 


"8 


5 


-0.178 


-7.78 


+12.8 


+43.7 


+103.8 


+67.9 


10 


-0.151 


-6.58 


+11.3 


+43.7 


+ 105.0 


+68.1 


15 


-0.145 


-6.33 


+ 11.1 


+43.7 


+105.5 


+68.3 


20 


-0.142 


-6.22 


+ 11.0 


+43.7 


+105.6 


+68.3 


25 


-0.142 


-6.20 


+ 10.9 


+43.7 


+105.6 


+68.4 



Table 5: Coefficients of the matrix elements that describe the short distance effects in a 
for ep = 0.012. 



cross section. Large positive values of the ^Pq matrix element are capable of creating a 
positive a. 

A notable feature which was observed in our numerical analysis is the dependence of 
a on Ep, which is very pronounced; |a| decreases with increasing sp. This corresponds 
to a similar behavior of a with the b quark mass as observed in ||3l|; the smaller (or 
larger ep) becomes, the smaller |a| gets to be. However, at the Tevatron, since the shape 
of the distribution function itself is not too important, this behavior is mainly related to 
determining the effective b quark mass, as introduced in (|5^). 

As stated earlier, in our analysis, the various matrix elements and ep are not com- 
pletely independent, but their values are correlated because of the conditions imposed in 
section |6.3| . Hence, the extreme values for a are a compromise of the tendencies described 
above. The maximum of a corresponds to values of {0(^^ {"^Si)) and {Og^^{^Si)) which are 
close to their lower limits as specified in appendix |A.1| . {Og^^{^So)) is also rather small, 
{Og^^{^Po)) is large and positive and ep is at its allowed maximum. The minimum of a is 
reached for Of^^l^Si)) and {Og^^ {"^Si)) being close to their upper limits, {Og^^{^So)) at 
an intermediate value and {O^^^ {^Pq)) being close to zero, but negative, ep is at its lower 
bound, thus illustrating the strong dependence of a on this parameter. 



7.3 Other theoretical uncertainties 

The variation of the factorization scale, /i, between its standard value mj, = 4.7 GeV and 
2mc = 3 GeV changes the results significantly. The obvious reason for this behavior is 
that lowering yU from 4.7 GeV to 3 GeV is equivalent to reducing the color singlet matrix 
elements by ~ 58% and at the same time increasing all octet matrix elements by ~ 6% as 
can be estimated from figure |1|. The net effect is a 10% reduction of the unpolarized cross 
section, da/dhx, and an increase of 20 . . . 25% in a. 

Increasing the light quark mass m/, g/ = d, s, from to 150 MeV causes a reduc- 
tion of the unpolarized cross section of approximately 1.5% if all other parameters are 
left unchanged. Again a is less strongly affected and increases by less than 1% . The 
uncertainties of the remaining quantities exclusively affect the cross section, because they 
have only an influence on its normalization and thus do not affect the predictions for a. 
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To this category, belong the CKM matrix element \Vcb\, the average B life time tb and 
the branching fraction Br{ip Ai+yU^) which is used to normalize the spectrum. In the 
combined error of da/dkT, these uncertainties have been included. 




5 10 15 20 25 30 

kj in GeV 



Figure 6: Unpolarized cross section da/dhx for J/tp production with contribution from all 
different channels and the combined error. 



7.4 Final results 

The final results for J/ ip as well as the experimental data are presented in figures ^ and |^. 
The corresponding quantities for ip', to which no feed-down channels contribute, are shown 
in figures ^ and p. 

The unpolarized J/ip cross section shows a good agreement with the experimental 
data 1^ within error bars at intermediate transverse momenta fc^. On the other hand it 



overestimates the data quite strongly for small and large /ct, a feature which is a result of 



extrapolating the B cross section (see section for details). 

Our prediction for the polarization parameter a is consistent with 0, but the central 
value prefers to be small and negative. For the higher transverse momentum range = 
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Figure 7: a for J/ip from B decays with contributions from all different channels and the 
combined error. 



10 ... 30 GeV, we have a = —0.04 ± 0.06. We find a good agreement with experimental 
data |p3, although one has to admit that the statistical experimental uncertainties are 



enormous and for some fcr bins cover as much as one third of the theoretically allowed 
parameter space of a. The prediction for a is almost independent of fcy and only shows a 
slight decrease towards the lowest transverse momentum values displayed in figure |^. Also 
the central values of the experimental data exhibit such a tendency, but the decrease at 
low kT is much more significant. Even here the predictions lie within la deviation of the 
experimental data. 

Obviously on the basis of the current data it is not possible to constrain the values of 
the poorly determined color octet ^Sq and ^Pq matrix elements any further. Nevertheless, 
there are good prospects that the situation on the experimental side will improve in the 
near future. Run II of the Tevatron is in progress for one and a half years already and the B 
physics group of the CDF collaboration expects to increase statistics effectively by a factor 
of 50 ||38|. This will reduce the statistical error at least in the low and medium transverse 



momentum range significantly and also data in higher kx bins will become available. Hence 
it is expected that after Run II only systematical errors will dominate the uncertainties 
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Figure 8: Unpolarized cross section da/dhx for ip' production with the combined error. 



of the a measurement Even if those will not be improved, a will have an error of 



the order of ~ 0.02, enough to exclude a good part of the numerical ranges of the matrix 
elements. 

Since our calculation can be equally applied to ifj' production without significant mod- 
ifications, we have extended the analysis to this charmonium state. For intermediate and 
large values of the unpolarized cross section in figure § agrees with the CDF data within 
error bars. The reason for the excess at low kx is the same as in the case of J / ip discussed 
earlier. With a = —0.03 ± 0.07 for kx = 10. . .30 GeV, the polarization parameter for 
ip' does not differ significantly from the one for J/tp, but here a comparison with data 
is almost impossible as can be seen from figure The situation is worse than for J/ip, 
because the error bars are huge and one of the data points even lies outside the allowed 
region for a, which is restricted to the interval [—1,+!] for theoretical reasons. On the 
other hand, once more precise data will be available, to derive tighter constraints on the 
ip' matrix elements will be more straight forward than for J/ip, because all ip' are directly 
produced from B decays, since no feed-down channels are known to contribute. 
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Figure 9: a for from B decays with the combined error. 

8 Summary and conclusions 

In this paper, using the NRQCD formahsm and the PM, we have calculated the semi- 
inclusive decay rate, V{B J/ip + X) with polarized J/ip as final states. Subsequently 
these results were generalized to the case of other J = 1 quarkonium states, i.e., ip' and 
Xci- The results were applied to the Fermilab Tevatron setup. In this case, we calculated 
the differential cross sections for unpolarized J /ip and ip' production in B decays and, the 
polarization parameter a for J /ip and ip' , originating from B decays. The B meson produc- 
tion cross section in pp collisions at the Tevatron was implemented by a phenomenological 



fit to the CDF data |3^. We considered the feed-down from ip' and Xci for J /ip production. 
To obtain a meaningful comparison of our predictions with experimental data, we carried 
out a detailed analysis of the various theoretical uncertainties involved. In particular, it 
was shown that a is almost not influenced by most input parameters, except for the ma- 
trix elements and the distribution function parameter, Sp. Therefore, for an extraction of 
this parameter, it is pertinent for more precise numerical values of the matrix elements, 
especially, the poorly determined (Og (^^o)) and (Og (^Pq))- 
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To perform such an improved fit of these non-perturbative matrix elements, the preci- 
sion of the data has to be increased significantly which is also expected in Run II of the 
Tevatron. Furthermore, it would also be desirable to experimentally separate direct J/ if) 
production from the feed-down channels. This restricts the uncertainties, simply because 
the number of relevant matrix elements get reduced . Theoretically, an inclusion of higher 
order corrections would be preferable to reduce the errors due to factorization scale depen- 
dence, in particular that of the color singlet Wilson coefficient. Since it is known that this 
cannot be achieved at next-to-leading order, a next-to-next-to-leading order calculation 
might be necessary Also, a better knowledge of would improve the precision of a. 



For if)' this could be accomplished with a fit to the unpolarized momentum spectrum of ip' 
in B decays to more accurate CLEO |72| and BaBar data that have become available 



very recently. The feed-down channels and resonant two body final states (J/V' + K and 
J/ip + K*) complicate such a fit for B J/ip + X. Additionally, it might be worth to 
try a different parameterization for the heavy quark distribution function, because in b 
fragmentation, which serves as a motivation for the distribution function, the Peterson 
form might not be appropriate . 



A significant reduction of the matrix element errors is mostly likely to be achieved with 
the help of a global fit. Among the processes that could contribute to this fit, a from B 
decays might play an important role, being one of the few quantities that is sensitive to 
the individual values of the {Og{^So)) and (Og (^Pq)) niatrix elements. 
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A Appendix 

A.l NRQCD Matrix Elements 

Here, we present numerical values of the matrix elements from the literature that are 
required in our analysis. The tables below are structured as follows. Apart from the 
numerical value, including statistical as well as systematical errors (wherever given) we 
refer to the method/process that has been used to extract them and the corresponding 
references to the original publications. This list only provides an overview of the numerical 
values that we have considered in our analysis and is by no means meant to be exhaustive. 
In the last line of each table we give a range for the corresponding matrix element that we 
have used in our calculation as described in section 
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A. 1.1 J/ip Matrix Elements 



{Oi'^^i^Si)) in GeV3 


metliod/process 


Reference 


0.994 ±0.0021°:^°^ 
1.16 

1.1 ±0.1 
0.763 ±0.054 
1.3±0.1 


lattice calculation 
Buclimiiller-Tye potential 

J/ip ^ e~^e~ 

J/tp e+e~ (without QCD) 
J/^p e+e- (incl. LO QCD) 




66 
74 

75 
75 




1.1 ±0.1 







(08^/'^(3^i)) in 10-3 GeV3 


method/process 


Reference 


6.6 ±2.1 
3.94 ±0.63 
10.6 ± 1.4+^°/ 
9.6 ±1.5 
4.4 ±0.7 


hadroproduction (CDF data) 
hadropro duct ion (CDF data) 
hadroproduction (CDF data) 
hadroproduction (CDF data) 
hadroproduction (CDF data) 


1 
1 
1 
1 
1 


67 
75 
36 
76 
18 


,38| 

1 

1 

1 

1 


6.0 ±3.0 







M//'^ in 10-2 GeV3 



method/process 



Reference 



6.6 ± 1.5 
6.52 ±0.67 
3.0 

8.7 ±0.9 
2.0 

4.38 ± 1.15: 
1.5 

1.32 ±0.21 



1.52 
■0.74 



3 

3.47 
7 

3.4 
7 

3.5 
3.1 

3 



hadroproduction (CDF data) 
hadroproduction (CDF data) 
hadroproduction (fixed target) 
hadroproduction (CDF data) 
photoproduction 
hadroproduction (CDF data) 
Br{B J/ij + X) (CLEO data) 
hadroproduction (CDF data) 



m 



m 



m 
m 

76| 



5.0 ±4.0 



3.5 



{O'J'^i'So)) in 10-2 GeV3 


method/process 


Reference 


4.0 
14.5 


leptoproduction 

Br{B J/i) + X) (CLEO data) 


1 
1 


20 
26 


1 
1 


6.5 ±5.5 







{Ol'^^i^P^)) in 10-2 GeV^ 


method / process 


Reference 


-0.3 ■ ml 
-5.51 


leptoproduction 

Br{B J/i) + X) (CLEO data) 


1 
1 


20 
26 


1 
1 


-4.0 ±6.0 
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A. 1.2 'ip' Matrix Elements 



(of in GeV^ 


method/process 


Reference 


0.440 ± 0.043 
0.65 ±0.6 
0.76 
0.11 


e+e- (incl. LO QCD) 
^' e+e" (incl. NLO QCD) 
Buchmiiller-Tye potential 
leptonic decay rate 


1 

1 
1 


26 
18 
74 
77 


1 

1 
1 


0.6 ±0.2 








{Ot\^Si)) in 10-3 GeV^ 


method/process 


Reference 


6.20 ±0.95 
4.2 

4.6 ± 1.0 
4.4±0.8l|^ 
4.2 ± 1.0 


hadroproduction (CDF data) 
hadroproduction (CDF data) 
hadroproduction (CDF data) 
hadroproduction (CDF data) 
hadroproduction (CDF data) 


1 
1 
1 

1 


26 
77 
67 

W\ 

18 


1 
1 

M 
i 


5.5 ±2.5 







Mf in 10-2 GeV=^ 



method / process 



Reference 



1.8 ±0.6 
1.79 ±0.51 
1.8±0.56+°j^ 
0.52 

1.3 ±0.5 
0.6 



3 

2.56 

3.5 

7 

3.5 
3.1 



hadroproduction (CDF data) 
hadroproduction (CDF data) 
hadropro duct ion 
hadroproduction (fixed target) 
hadroproduction (CDF data) 
Br{B + X) (CLEO data) 



m 



m 







1.75 ± 1.25 



3.0 



(Cg^'V^o)) in 10-2 GeV^ 


method/process 


Reference 


-0.96 


Br{B J/i) + X) (CLEO data) 




26 


1 


2.5 ±2.5 








(Cf (3Po)) in 10-2 GeV^ 


method/process 


Reference 


2.58 


Br{B J/i) + X) (CLEO data) 




26 


1 


0.0 ±10.0 
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A. 1.3 Xcj Matrix Elements 



(Cf°(3Po)) in 10-2 GeV^ 


method/process 


Reference 


4.8 ■ ml 
4.4 ■ ml 
22.9 ±2.5 

8.8 ±2.13 

8.9 ±1.3 


Buclimiiller-Tye potential 
Buchmiiller-Tye potential 
hadropro duct ion (CDF data) 
hadronic Xcj decays 

Xc2 77 


1 

1 
1 


74 
74 
75 
25 
18 


44| 


6.0 ±4.0 







in 10-3 GeV^ 


method / process 


Reference 


4.5... 6.5 
0.681 ±0.175 

3.2 ± 1.4 
1.39 ±0.17 

2.3 ±0.3 


B^Xc2 + X 

hadroproduction (CDF data) 
r{B^Xc2 + X)/r{B^eu, + X) 
hadroproduction (CDF data) 
hadroproduction (CDF data) 




41 
75 
77 
26 
18 


1 
1 
1 


3.5 ±3.0 
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